







%Clculo de la masa de la Volante por el mtodo de Wittenbauer
%Datos de entrada


clc,clear all, close all
l1=0.05; l2=0.25; lbs2=0.10;
D1=0.13; D2=0.11;
J0=0.07;
J2=0.025; m2=1.8; m3=2.2;
omegamed=(pi*800)/30;
deltaprim=0.0125;

Ndatos=3600;

fi1=linspace (0,360,Ndatos+1)
fi1=fi1/180*pi
carrera=l1*2;
pms=l1+l2
pmi=l2-l1















%Cinemtica del mecanismo

omega1=1; %slo para la reduccin
for N=1:Ndatos+1
xc(N)=l1*cos(fi1(N))+(l2^2-l1^2*sin(fi1(N))^2)^(1/2);
sc(N)=(xc(N)-pmi)/carrera;
fi2(N)=-asin(l1*sin(fi1(N))/l2);
omega2(N)=-l1*cos(fi1(N))*omega1/l2/cos(fi2(N));
vxc(N)=-l1*omega1*sin(fi1(N)-fi2(N))/cos(fi2(N));

vxb(N)=-l1*sin(fi1(N))*omega1;
vyb(N)= l1*cos(fi1(N))*omega1;
vxs2(N)=vxb(N)-lbs2*sin(fi2(N))*omega2(N);
vys2(N)=vyb(N)+lbs2*cos(fi2(N))*omega2(N);
vs2(N)=sqrt(vxs2(N)^2+vys2(N)^2);
end








%Diagramas indicados

for N=1:Ndatos+1 
    if and(sc(N)<=1.0 && sc(N)>=0.87  , fi1(N)>=0 && fi1(N)<=pi)==1
         pI(N)=0.7226e7*sc(N)^2-1.1782e7*sc(N)+0.4781e7;       
    elseif and(sc(N)<0.87 && sc(N)>0  , fi1(N)>0 && fi1(N)<pi)==1
        pI(N)=0;
    elseif and(sc(N)>=0 && sc(N)<=0.627 , fi1(N)>=pi && fi1(N)<=2*pi)==1
        pI(N)=1.6940e5*sc(N)^2+2.5264e5*sc(N);
    else
         pI(N)=2.25e5;
    end   
end

plot(sc,pI),title('Diagrama Indicado del compresor'),xlabel ('s'),ylabel 				('presin Pa')

 for N=1:Ndatos+1 
     if and(sc(N)<=1.0 && sc(N)>=0.391   , fi1(N)>=0 && fi1(N)<=pi)==1
            pII(N)=0.3955e6*sc(N)^2-1.2890e6*sc(N)+1.1185e6;
      elseif and(sc(N)<0.391 && sc(N)>0  , fi1(N)>0 && fi1(N)<pi)==1
            pII(N)=6.75e5;
      elseif and(sc(N)>=0 && sc(N)<=0.258 , fi1(N)>=pi && fi1(N)<=2*pi)==1
             pII(N)=2.7547e6*sc(N)^2-2.4549e6*sc(N)+0.6750e6;              
      else 
         pII(N)=2.25e5;
     end   
 end

hold on
 plot(sc,pII);grid






%Fuerzas sobre el cilindro
FI=pi*D1^2/4*pI;
FII=pi*(D1^2-D2^2)/4*pII;
Fp=FII-FI;










%Reduccin del Momento
for N=1:Ndatos+1
    cosenoalfa=(sign(Fp(N))*sign(vxc(N)));
    Mc(N)=-abs(Fp(N))*abs(vxc(N))*cosenoalfa;%ojo con el "menos" y con 									los valores absolutos 
end
Mc

figure
plot(fi1,Fp),title('Fuerza en el cilindro'),xlabel ('fi1 rad'),ylabel 										('fuerza N')
grid
figure
plot(fi1,Mc),title('Momento reducido - Momento motor'),xlabel ('fi1 						rad'),ylabel ('Momento N.m')
grid









%Clculo del Momento Motor
%Integracin total del momento por trapecios 

deltafi=fi1(2)-fi1(1);
Area=0;
for N=1:Ndatos
   A(N)=(deltafi*Mc(N))+0.5*(deltafi*(Mc(N+1)-Mc(N)));
   Area=Area+A(N);
end

Mm=Area/(2*pi)
Mmot=ones(1,Ndatos+1)*Mm;
hold on
plot (fi1,Mmot,'r')








%Clculo de la Momento variacin de la energa cintica
%Integracin por trapecios 

EK(1)=0;
for N=2:Ndatos
   DeltaEK(N)=(deltafi*Mmot(N))-deltafi*(Mc(N)+0.5*(Mc(N+1)-Mc(N)));
   EK(N)=EK(N-1)+DeltaEK(N); 
end
EK(Ndatos+1)=EK(1);
EK;
figure
plot(fi1,EK),title('Incrementos de Energa Cintica'),xlabel ('fi1 						rad'),ylabel ('Energa J')
grid










%Reduccin del Momento de inercia de los eslabones sin incluir manivela
for N=1:Ndatos+1
   Jeslab(N)=J2*omega2(N)^2+m2*vs2(N)^2+m3*vxc(N)^2;
end
figure
plot(fi1,Jeslab),title('Momento de inercia reducido'),xlabel ('fi1 						rad'),ylabel ('Momento Inercia kg.m^2')
grid







%Diagrama de Wittenbauer
figure
plot(Jeslab,EK),title('Diagrama de Wittenbauer'),xlabel ('Momento Inercia 				kg.m^2'),ylabel ('Energa J')
grid










%%%%%%%%%%%%%%% determinacion de las tangentes max min

format long
tanksimaxprim=(1/2*omegamed^2*(1+deltaprim))
tanksiminprim=(1/2*omegamed^2*(1-deltaprim))

ksimaxprimgrad=(atan(tanksimaxprim))*180/pi
ksiminprimgrad=(atan(tanksiminprim))*180/pi
format

p=[tanksimaxprim 20]; %recta de pendiente "tanksimaxprim" y desplazada en 				20 unidades
xksimax =[2e-3 11e-3];
yksimax =polyval (p,xksimax);
hold on
plot (xksimax,yksimax,'k')

dmin=1e6
 for N=1:Ndatos+1
     P=[Jeslab(N) EK(N)];   
    Q1= [xksimax(1) yksimax(1)];
    Q2= [xksimax(2) yksimax(2)];
    d(N)= abs(det([Q2-Q1;P-Q1]))/norm(Q2-Q1);
    if d(N)<dmin
      dmin=d(N);
      Nmin=N;
    end
 end
 
b2=-tanksimaxprim*Jeslab(Nmin)+EK(Nmin)
p2=[tanksimaxprim b2]; %recta tangente a Wittenbauer con pendiente 							tanksimaxprim

yksimax =polyval (p2,xksimax);

plot (xksimax,yksimax,'r')
plot ( Jeslab(Nmin),EK(Nmin),'*r')








%%%%%%%%%%%%%%%%%%
p=[tanksiminprim -80]; %recta de pendiente "tanksiminprim" y desplazada 						en -80 unidades
xksimin =[2e-3 11e-3];
yksimin =polyval (p,xksimax);
plot (xksimin,yksimin,'k')

dmin=1e6
 for N=1:Ndatos+1
     P=[Jeslab(N) EK(N)];   
    Q1= [xksimin(1) yksimin(1)];
    Q2= [xksimin(2) yksimin(2)];
    d(N)= abs(det([Q2-Q1;P-Q1]))/norm(Q2-Q1);
    if d(N)<dmin
      dmin=d(N);
      Nmin=N;
    end
 end
 
b3=-tanksiminprim*Jeslab(Nmin)+EK(Nmin)
p3=[tanksiminprim b3]; %recta tangente a Wittenbauer con pendiente 											tanksiminprim

yksimin =polyval (p3,xksimin);

plot (xksimin,yksimin,'r')
plot (Jeslab(Nmin),EK(Nmin),'*r')











%solucin del sistema de ecuaciones lineales
%p2=[tanksimaxprim b2];
%p3=[tanksiminprim b3];

format long
A=[tanksimaxprim -1;tanksiminprim -1]
B=[-b2;-b3]
c=A\B
format

%Magnitud del momento de inercia de la volante
Jvol=-c(1)-J0
1


